Predicting Ion Diffusion from the Shape of Potential Energy Landscapes

We present an efficient method to compute diffusion coefficients of multiparticle systems with strong interactions directly from the geometry and topology of the potential energy field of the migrating particles. The approach is tested on Li-ion diffusion in crystalline inorganic solids, predicting Li-ion diffusion coefficients within 1 order of magnitude of molecular dynamics simulations at the same level of theory while being several orders of magnitude faster. The speed and transferability of our workflow make it well-suited for extensive and efficient screening studies of crystalline solid-state ion conductor candidates and promise to serve as a platform for diffusion prediction even up to the density functional level of theory.


Choice of single-particle grid sizes
The multi-particle MMC routine in the Ionic TuTraSt workflow was run both with a 0.1 Å × 0.1 Å × 0.1 Å and 0.2 Å × 0.2 Å × 0.2 Å single-particle grid as input to test the effect of the grid size resolution on the stability and accuracy.Figure S1 shows that the accuracy of the the prediction of diffusion coefficients compared with MD does not decrease when the single-particle grid resolution decreases from 0.1 Å to 0.2 Å.Rather, the finer grid resolution shows a lower stability with three false positives as the Li-counts are spread over a larger number of volumetric bins, thus rendering more noise in the statistics and a less smooth PES surface compared to the 0.2 Å grids.

Effect of Li partial charge for validation set 2
Different partial charges were tested for Li during the multi-particle MMC simulations; q = 0.25e, 0.5e, 0.75e and 1e, as well as Li charges computed by the REPEAT method as   In the third case, for the LiNiN structure, the false positive result is observed from the Arrhenius plot (Figure 6) at lower temperatures.At 1000 K the agreement is accurate, however from 700 K and below, a long range ordering in the x-direction emerges that extends across the full simulation.This long-range ordering breaks the symmetry and the basins do not merge into percolated diffusion channels.This effect is not seen in the multi-particle PES as the MMC simulation does not get stuck in local minima and can sample all low energy states more efficiently.Also, when constructing the multi-particle PES the Li sampling counts are averaged over all unit cells to a single unit cell PES, and thus possible long-range effects are not detected.In contrast, at 700 K, the periodic unit in the x-direction extends across the whole simulation cell, consisting of two inverted blocks with four unit cells each (shaded and transparent boxes, respectively).

Efficiency dependencies
Figure S7 shows how the speed-up ratios (as defined by the order analysis), of the respective single-particle grid and multi-particle MMC routines vs MD sampling, depend on the stoichometric Li-ratio and the number of grid points.Here, the dependencies are presented for different scaling orders; O(N 2 ) for direct summation of the Ewald method, O(N 3/2 ) for standard Ewald summation as implemented in LAMMPS (used in this study) and O(N log N ) corresponding to the particle-particle particle-mesh (PPPM) method.S6 When using classical force-fields for the both single-particle grid sampling as well as the multi-particle MMC routine, the computational cost of the MMC routine dominates over the single-particle grid construction.The order of speed-up of Ionic TuTraSt over

Figure S1 :
Figure S1: Directional diffusion coefficients for validation set 1 computed with singleparticle grid sizes 0.2 Å (left) 0.1 Å (right) on the y-axis relative to the corresponding diffusion coefficients computed with MD on the x-axis on a log-log scale.The dashed lines guide the limits for deviation of one and two orders of magnitude, respectively.

Figure S2 :
FigureS2: Comparison of Ionic TuTraSt performance using different charges (q) for the Li during the multi-particle MMC simulation for validation set 2. Plots on the left show diffusion coefficients obtained from Ionic TuTraSt against from diffusion coefficients obtained from PBMD simulations.The dashed lines guide the limits for deviation of one and two orders of magnitude, respectively.Histograms on the right show the difference in break-through (BT) energies obtained by Ionic TuTraSt relative to to PBMD.

Figure S4 :
FigureS4: Potential energy isosurfaces in the xy-plane at 20 kJ/mol from MD simulations for structure Li 3 YBi 2 (in white).The configurations that the three stoichiometric Li ions can take in the unit cell are presented in red.The radial distribution function of the MD has a first peak at around 2.5 Å, which is the approximate lengths between the sides occupied in the shown configurations (red), suggesting that it is highly improbable that sides with lower distance between them are simultaneously occupied.

Figure S5 :
Figure S5: The Li−Li radial distribution function in Li 3 YBi 2 computed from MD simulations.

Figure S6 :
Figure S6:The MD PES isosurface (20 kJ/mol) of LiNiN from the 700 K simulation (right) shows a long-range Li + ordering compared to 1000 K (left).The dashed lines show the 8 × 8 unit cell boxes building up the simulation cell in the xy plane.The red boxes show the repeating unit for the respective temperature.At 1000 K, the repeating unit corresponds to the unit cell.In contrast, at 700 K, the periodic unit in the x-direction extends across the whole simulation cell, consisting of two inverted blocks with four unit cells each (shaded and transparent boxes, respectively).
FigureS7: Speed-up dependencies presented for different scaling orders.Left: Log-scale speed-up ratios of MMC compared to MD depending on the Li-ratio.Right: Log-scale speed-up ratios of the grid sampling routine compared to MD depending on the number of 0.2 Å grid points.